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1.0  INTRODUCTION 


A  modem  nuclear  submarine  radiates  acoustic 
energy  as  it  moves  underwater  in  the  ocean.  It  is  possible  to 
detect  and  track  these  submarines  by  using  this  radiated 
energy.  Typically,  this  energy  is  acquired  by  passive  sonobuoys 
and  processed  electronically  to  give  estimates  of  the  Doppler- 
shifted  frequency  associated  with  the  submarine  and  also  an 
estimated  bearing. 

1 . 1  Background  -  Processing  Alternatives 

Once  the  raw  data  from  the  sonobuoys  has  been 
processed  to  give  frequency  and  bearing  estimates,  it  is  passed 
through  a  tracking  algorithm  which  will  produce  estimates  of 
the  target’s  state  vector,  i.e.,  position,  velocity,  etc. 

At  present  there  are  two  primary  classes  of 
algorithms  for  solving  this  problem- -batch  algorithms  and 
sequential  algorithms.  The  names  are  derived  from  the  way  each 
algorithm  processes  data.  Batch  algorithms  simultaneously 
process  several  data  points  (i.e.,  a  batch  of  data)  at  once  to 
provide  tracking  parameters,  while  sequential  algorithms  process 
one  data  point  at  a  time,  giving  updated  state  vector  estimates 
after  each  data  point  has  been  processed.  Over  the  past  several 
years,  Tracor  has  developed  several  batch  algorithms  called 
Maximum  Likelihood  Procedures  (MLP)  to  track  submarine  targets 
for  a  number  of  specialized  tracking  system  configurations. 

Each  of  these  procedures  has  been  tested  on  sets  of  real-world 
and  simulated  data  with  varying  degrees  of  success.  However, 
in  doing  so  several  points  became  clear: 
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(1)  The  MLP  is  self -initializing .  That  is, 
good  initial  estimates  of  the  state  vector 
can  be  generated  from  the  data.  There  is 
no  need  to  supply  an  accurate  initial  guess. 

(2)  The  MLP  requires  several  different  motion 
models  to  characterize  possible  submarine 
maneuvers.  Model  selection  can  become  a 
problem  in  terms  of  time  and  appropriateness. 

(3)  Because  the  MLP  processes  data  in  a  batch, 
extremely  accurate  estimates  of  the  state 
vector  can  be  obtained  when  the  proper 
motion  model  has  been  selected.  However, 
serious  loss  of  track  can  occur  when  the 
model  currently  being  used  no  longer 
becomes  valid  but  the  selection  process 
has  not  yet  selected  a  new  model. 

(4)  The  algorithm  possesses  moderate  to 
substantial  computer  time  requirements, 
depending  upon  the  complexity  of  the  track 
being  attempted,  and  moderate  storage 
requirements . 

To  provide  a  means  of  comparison  for  the  MLP 
and  also  in  an  attempt  to  gain  insight  into  the  nature  of 
sequential  algorithms,  Tracor,  in  conjunction  with  Dr.  Byron 
Tapley  of  the  University  of  Texas,  designed  and  implemented 
a  sequential  tracking  algorithm  based  on  the  extended  Kalman 
filter.  While  testing  this  algorithm  several  things  became 
clear : 
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(1)  This  algorithm  has  smaller  core  storage 
and  execution  time  requirements  than  the 
MLP . 

(2)  Tracking  accuracy  is  as  good  or  better 
than  the  MLP. 

(3)  The  sequential  algorithm  is  highly 
susceptible  to  initialization  errors, 
requiring  a  fairly  close  initial  state 
vector  estimate  to  maintain  good  track. 

(4)  The  sequential  model  can  be  implemented 
with  a  single  stochastic  motion  model. 

This  eliminates  the  need  for  the  model 
selection  process,  and  also  allows  some 
"slack"  in  estimating  the  tracking 
parameters . 

1.2  Hybrid  Research  Effort  and  Program  Objectives 

After  testing  of  the  MLP  and  sequential 
algorithms  revealed  their  complementary  strengths,  it  was  then 
decided  that  a  combination  of  the  two  approaches,  or  a  hybrid 
algorithm,  should  be  developed  to  take  advantage  of  these 
strengths.  It  was  hoped  that  the  MLP  algorithm  could  intialize 
the  tracking  procedure  for  the  sequential  algorithm  and  also 
provide  reinitialization  should  the  sequential  filter  lose 
track.  A  preliminary  version  of  the  hybrid  was  created  using 
the  MLP  and  sequential  algorithms  and  early  tests  showed 
potential.  However,  several  problems  remained.  Among  them 
were : 
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(1)  The  MLP  was  too  cumbersome  to  use 
efficiently  as  an  initializer;  a  more 
efficient  batch  procedure  was  needed. 

(2)  Information  transfer  between  the  batch  and 
sequential  portions  had  to  be  streamlined 
and  made  more  efficient. 

(3)  The  sequential  algorithm's  numerical 
properties  needed  to  be  improved  for 
smaller  word  length  computers. 

(4)  The  capability  for  handling  maneuvering 
sensors  was  needed. 

(5)  The  capability  for  tracking  multiple 
targets  was  needed. 

(6)  Improvement  and  refinement  of  the  switching 
rules  from  the  batch-to-sequential  and 
sequential -to -batch  modes  of  operation 
were  needed. 

1 . 3  Summary  of  Results 

A  hybrid  algorithm  incorporating  the  above 
features  has  been  designed  and  implemented.  This  study 
provides  a  comparison  between  the  hybrid  algorithm  and  the  MLP 
and,  in  addition,  a  comparison  with  a  sequential  algorithm 
that  contains  an  iterated  sequential  starter. 
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THEORY  AND  EXPLANATION  OF  ALGORITHMS 


Maximum  Likelihood  Procedure  (MLP) 


2.1.1 


Motion  Models 


The  MLP  algorithm  is  based  on  the  assumption  that 
submarines  typically  execute  a  fairly  small  set  of  maneuvers, 
including : 

(1)  Constant  Linear  Velocity 

(2)  Constant  Linear  Acceleration 

(3)  Constant  Radius  and  Angular  Velocity  Turns 

(4)  Constant  Radius  and  Angular  Acceleration  Turns 

f 

Each  of  these  maneuvers  has  a  corresponding 
motion  model.  They  are: 

(1)  Constant  Velocity 


X(t)  ■  X„  +  V  t 
o  x 


Y(t)  -  Yq  +  Vyt 


V*>  -  Vx 


vyct)  -  vy 
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(2)  Constant  Acceleration 

X(t)  =  X0  +  Vxt  +  .5  axt> 

Y(t)  =  Yq  +  Vyt  +  .5  ayt2 

"  Vx  +  axC 
Vy(t)  -  Vy  +  ayt 

(3)  Constant  Speed  Turn 

X(t)  *  X_  +  V  sin(wt)/w  -  V  [1  -  cos(vt)]/w 
ox  y 

Y(t)  =  Y  +  V  (1  -  cos(wt)]/w  +  V  sin(wt)/w 
u  x  y 

vx(t)  -  Vxcos(wt)  -  Vysin(wt) 

V  (t)  =  Vxsin(wt)  +  Vycos(wt) 

(4)  Constant  Angular  Acceleration,  Constant 
Radius  Turn 

X(t)  =  XQ  +  {Vxsinwt (1+. 5at)  -  V  [ 1-coswt (1+. 5at) ] } /w 
Y(t)  =  Y  +  {V  [1-coswt  (l+.5at)  +  V  sinwt  (1+.  5at)  ] }  /w 

v  a  y 

Vx(t)  =  (1+at) (Vxcoswt (1+. 5at)  -  Vysinwt (1+. 5at) } 

V  (t)  =  (1+at) {Vxsinwt (1+. 5at)  +  V„coswt (1+. 5at) } 

y  x  y 
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Each  of  these  maneuvers  is  characterized  by  a 
unique  set  of  target  parameters.  For  example,  constant  velocity 
motion  is  characterized  by  the  following  vector: 


where 


Y  = 


Xt  =  X-position  at  time,  t 

Yt  =  Y-position  at  time,  t 

Vv  =  X-velocity  at  time,  t 
t 

Vy  =  Y-velocity  at  time,  t 

*,fn  =  Source  frequency  parameter  for  each 
target  source  frequency 


The  approach  used  by  MLP  to  track  targets  is 
summarized  in  Figure  2-1.  Frequency  and  bearing  data  from  each 
sonobuoy  maintaining  contact  with  the  target  are  input  to  each 
maneuver  model  for  which  there  is  a  non -zero  a  priori 
probability  that  the  target  could  be  executing  that  maneuver. 

The  Maximum  Likelihood  procedure  is  applied  to  estimating  the 
parameters  associated  with  each  model;  these  parameter  values 
and  the  residual  data  errors  are  compared  for  each  of  the  models 
and  three  decisions  are  reached. 


(1)  The  maneuver  model  and  parameter  vector 
associated  with  that  model  are  selected. 
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FIG.  2.1 


PROGRAM  FLOW  FOR  MANEUVERING  TARGET  TRACKING 
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(2)  The  maneuver  models  which  are  totally 
inadequate  to  explain  che  data  are 
determined. 

(3)  The  reduced  tracking  intervals  for  each 

of  the  models  which  could  explain  the  data 
are  selected. 

At  the  conclusion  of  this  procedure  new  data 
are  input  and  the  procedure  is  repeated.  Note  that  the  MLP 
passes  the  data  through  four  Maximum  Likelihood  parameter 
estimation  algorithms ,  each  one  corresponding  to  one  of  the 
maneuver  models  described  above. 

2.1.2  MLP  Algorithm 

The  mathematical  structure  of  the  Maximum 
Likelihood  procedure  is  the  same  regardless  of  the  motion  model 
and  data  types  used.  Associated  with  each  measurement  is  a 
nonlinear  function  that  generates  an  estimate  of  that  data 
measurement  given  the  time  of  the  measurement,  a  set  of  motion 
parameters,  and  a  motion  model.  For  example,  if  a  frequency 
measurement  is  taken  using  sensor,  j,  at  time,  t^,  the  estimated 
frequency  is  given  by, 

f(ti,X)  =  f0/[l+V(ti)-Rj(ti)/c|Rj(ti)|  ]  , 


where 


X  =  motion  parameter  vector, 

V(t.)  =  velocity  vector  of  target  at  time,  t^, 
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Rj(ti)  =  range  vector  of  the  target  from  sensor, 
j ,  at  time ,  t^ 

f0  *  center  frequency 
c  =  velocity  of  propagation  of  sound 

If  a  bearing  measurement  is  considered,  then 

f  (t  i ,  y)  =  tan'1  [Yj(ti)/Xj(ti)] 


where 


Xj(ti)  =  X-coordinate  of  target  at  time,  t^, 
relative  to  sensor,  j 

Yj (t^)  =  Y-coordinate  of  target  at  time,  t^, 
relative  to  sensor,  j 

The  program  can  also  use  range,  Doppler  ratio,  Doppler 
difference,  and  time  difference  of  arrival.  The  appropriate 
motion  model  is  used  to  generate  Rj (t)  and  Vj (t)  for  any  time, 
t. 


Let  Ej_  -  Zi  -  f(t^,y).  It  is  desired  that  the 
f's  be  close  to  the  measurement  values  Zj_.  A  reasonable  measure 
to  consider  for  this  closeness  is  the  sum  of  squares  of  residuals 
for  a  given  parameter  vector, 

S(y)  =  l  E ±2/a±2, 
i=l 

where  a^2  is  the  measurement  error  of  the  observation,  Zj_.  If 
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the  above  sum  is  small,  the  f's  will  be  close  to  the  Z^'s. 

If  it  is  assumed  that  the  measurement  values  are  perturbed  from 
the  true  mean  values  of  f(tj_,y)  by  Gaussian  random  variables, 
the  minimization  of  the  sum  is  equivalent  to  maximizing  the 
likelihood  function  L(Q)  of  the  sample  Z  as 


L(n)  =  EXP 

N 

-  1  Ei2/cJi2 

/ 

(2tt)N/2  t?  at 

i=l 

i-1 

Even  if  the  measurement  errors  are  not  Gaussianly  distributed, 
this  procedure  provides  a  least  squares  fit  to  the  data.  Thus, 
the  problem  statement  is:  Given  a  set  of  measurements,  Z,  and 
a  motion  model,  find  the  motion  parameter  vector,  y,  that 
minimizes 

N 

I  Ei 2 /ai2  . 

1=1 


Since  f(tj_,y)  is  a  nonlinear  function  of  the 
vector  ,y,  a  direct  solution  is  not,  in  general,  possible. 
However,  a  procedure  may  be  used  in  which  a  linear  approximation 
is  applied  iteratively  to  search  for  a  solution.  If  the  non¬ 
linear  problem  is  sufficiently  well  behaved  the  iterative 
technique  will  converge  to  the  desired  solution.  In  general, 
it  is  not  possible  or  is  too  difficult  to  prove  whether  this 
technique  will  converge  in  a  given  problem.  Hence,  the 
algorithm  is  executed.  If  it  converges,  a  solution  is  found; 
if  not,  add  more  data  and  try  again. 

To  linearize  the  problem  each  function,  f,  is 
expanded  in  a  first  order  Taylor  series  about  an  assumed 
parameter  vector,  y0 ,  as, 
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f(ti,Y)  =  f(ti,Y0)  + 


H 


l  ski 

k-i  3Yk 


•i  »Y0 


where  <5k  is  the  variation  of  the  k—  parameter  of  y.  If  the 
vector,  E0,  is  defined  by  the  elements, 

E0  -  zi  '  f(ti,Y0) 
and  the  matrix,  X,  by  the  elements 


3f(t  ,Y) 
3Yj 


t 


o 


Then,  the  linear  approximation  problem  is 

X5  =  Eg  . 


The  vector,  5,  may  be  found  in  the  usual  way  for  minimizing 
the  weighted  sum  of  squares  of  residuals  (Reference  1) .  The 
normal  equations  are 

XTa_1X<5  =  XTa_1  e0 


and  6  is  found  to  be 

5  =  (XTa"1X)-1XTa'1e0  . 


The  new  motion  parameter  vector  is  given  by 
Y  -  Y0  +  <$, 
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and  when  the  change  in  each  component  of  the  motion  parameter 
vector  becomes  small,  the  algorithm  is  deemed  to  have  converged 
and  iteration  stops.  The  covariance  of  the  estimate  y  is  given 
by  the  following  equation, 

P  =  (XTa'1X)"1  . 

One  of  the  key  features  of  the  MLP  is  the  method 
employed  to  invert  the  symmetric,  positive,  semi -definite  matrix, 
X^a-1X.  The  usual  inversion  method  is  to  use  a  Choleski 
decomposition,  that  is,  to  find  an  upper  triangular  matrix,  U, 
such  that  UTU  •  X^a  X.  The  matrix  U  is  easily  inverted  and 
U_1U~T  is  the  required  inverse  for  finding  the  change  in  y.  It 
would  be  desirable  to  find  U  directly  and  then  the  cross-product 
array  would  not  have  to  be  formed. 

Such  a  procedure  is  possible  if  one  uses  a 
sequence  of  Givens'  rotations  (Reference  2)  to  decompose  a  ^X 
into  the  product  of  an  orthogonal  matrix,  Q,  and  an  upper 
triangular  matrix,  U.  Then 

XTo-1X6  =  XTa-1e0 


becomes 


UTQTQU  =  UTQTa"^e0. 

Since  Q  is  orthogonal,  Q^Q  =  I  and  if  U  has  an  inverse,  then 
the  above  becomes 
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This  triangular  system  of  equations  is  easily  solved,  and  the 
covariance  matrix  is 

p  =  iTWt  . 

The  Givens'  notation  is  carried  out  on  each  row 
of  X;  hence,  it  is  only  necessary  to  store  U  and  one  row  of  X 
at  a  time.  Figure  2.2a  depicts  the  usual  Givens’  method  approach. 
Note  that  M  is  the  number  of  parameters  estimated  and  N  is  the 
number  of  data  points.  Because  this  procedure  is  computationally 
expensive,  there  is  an  alternative  procedure  developed  by 
Gentleman  (Reference  3)  which  is  used  in  the  MLP  algorithm.  It 
is  depicted  in  Figure  2.2b,  and  does  not  require  any  square  roots 
and  reduces  the  number  of  multiplications  by  one -quarter. 

Thus,  the  revised  algorithm  calls  for  the 
accumulation  of  D,  the  diagonal  matrix;  £,  an  upper  unit 
triangular  matrix;  and  6,  a  vector  of  transformed  residual 
errors.  The  array  of  weighted  derivatives,  o  ^X,  is  replaced 
by  QD2U,  and  a  2eQ  is  replaced  by  D  ^Q'a  which  is  equal  to 

0.  Using  the  normal  equations,  this  yields: 

XTa-1XS  =  XTa-1 eQ 

UTdVqD%  =  trTD^QTa“^G0 

DCS  =  D*QTa’*eQ 


U<5  =  0  . 
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CONSIDER  TWO  ROW  VECTORS 


0  .  .  . 

0  .  .  . 


0  \  ’W 


O  xt  x1+1 


Transform  by 

u‘k  -  cuk  +  sxk 

x'k  "  "S^k  +  Cxk 


where 


REQUIRING  ABOUT  2NM2  MULTIPLIES  AND  NM  SQUARE  ROOTS 


FIG.  2.2a  -  GIVENS'  ROTATIONS 
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CONSIDER  TJ  -  D^U 

D  *  Diagonal 

R  ■  Unit  Upper 

Triangular  Matrix 

ROTATE 

0  ...  0  /3  . . .  /3  Uk 

0  ...  0  /6x^  ...  xk 

Into 

0  ...  0  /S'  ...  /S  u’k 

0  ...  0  0  ...  /T’x  ’  k 

d  +  <5x? 

d6/d ' 
d/d' 

6x^/d' 

*k  "  xiUk 
Crk  +  Sxk 

REQUIRING  ABOUT  3/2  NM2  MULTIPLIES 

FIG.  2.2b  -  SQUARE  ROOT  FORMULATION  OF  GIVENS’  ROTATIONS 
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The  reduction  of  the  weighted  residual  errors  to  0  is 
accomplished  by  considering  them  as  an  augmenting  column  of  the 
matrix  X.  The  last  equation  involves  a  unit  upper  triangular 
matrix  which  makes  the  solution  for  6  particularly  easy. 

To  determine  which  motion  model  to  select,  the 
MLP  uses  a  statistical  test,  based  on  the  residual  errors 
generated  by  each  model,  called  an  F-test.  Under  the  normal 
distribution  of  errors  assumption,  the  residual  sums  of  squares 
from  each  model  are  distributed  as  chi-square  random  variables 
and  their  ratios  are  distributed  as  F-random  variables. 
Strictly,  the  F-test  may  only  be  applied  when  the  errors  are 
normally  distributed,  however,  it  has  been  shown  to  be  a  robust 
test  against  other  distributions  as  long  as  they  are  symmetric 
about  zero.  The  use  of  the  ratio  is  valuable  because  it  makes 
the  test  independent  of  scale  changes  in  the  input  variance  of 
the  data. 


In  selecting  a  model  it  is  desirable, 
intuitively,  to  favor  the  constant  velocity  model  over  the 
others  because : 

(1)  It  is  the  simplest  of  the  four  models. 

(2)  The  other  three  models  yield  good 
solutions  to  a  constant  velocity  data 
interval . 

This  is  done  by  conducting  a  hypothesis  test  at  a  level,  a. 
The  test  statistics  are: 
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Fij  =  (SS./DF.J/CSSj/DFj) 
i. j  =  1,2, 3, 4;  i<j 

where  SS^  and  DF^  are  the  stun  of  squares  of  residuals  and  the 
degrees  of  freedom  for  the  i—  model. 

The  hypothesis  test  is: 

If  Pr (f>F1^ )  <  a 

j  =  2,3,4;  accept  Model  1 


otherwise 


If  Pr(f>Fi_.)  <_  .5  for  all  j>i; 
i  =  2,3,4;  accept  Model  i 

It  has  been  found  that  an  a  of  about  .2  works  well.  This  value 
gives  somewhat  of  an  edge  to  Model  1  over  Models  2,  3,  and  4. 

The  choice  between  Models  2,  3,  and  4  is  equal  because  the 
models  are  of  roughly  equal  complexity. 

An  additional  statistical  test  is  performed  to 
determine  the  proper  data  interval  for  the  given  motion  model. 
Essentially,  the  test  examines  the  residual  stream  from  each 
data  source  for  a  mean  shift,  which  is  assumed  to  indicate  that 
the  target  has  undertaken  some  maneuver  which  cannot  be  described 
by  the  current  motion  model.  When  this  shift  is  detected,  a  new 
data  interval  is  begun  and  a  new  model  is  used. 
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The  above  exposition  is  not  intended  to  be  an 
exhaustive  examination  of  the  MLP  algorithm.  This  is  contained 
in  Reference  4.  Indeed,  the  MLP  is  significantly  more  complex 
and  possesses  far  more  detail  than  has  been  indicated  here. 

For  example ,  the  MLP : 

(1)  Uses  the  Marquardt  algorithm  and  several 
additional  rules  to  speed  convergence  to 
the  nonlinear  solution; 

(2)  Provides  and  uses  observability  measures 
and  accuracy  measures  for  solutions 
generated  by  the  algorithm; 

(3)  Provides  for  continuity  constraints  when 
switching  from  one  motion  model  to  another; 

(4)  Provides  for  the  input  of  a  priori  state 
vector  and  covariance  information; 

(5)  Provides  an  outlier  accommodation  scheme 
which  has  been  incorporated  into  the 
measurement  model . 

Rather,  the  intent  of  this  section  is  to  give  a  brief  overview 
of  the  main  features  of  the  MLP  and  to  explain  certain  ideas 
which  are  used  in  both  the  MLP  and  hybrid  algorithms. 

2 . 2  Hybrid  Algorithm 

2.2.1  Introduction 


The  current  Tracor  hybrid  algorithm  is  composed 
of  two  distinct  entities--a  scaled  down  MLP-type  algorithm  and 
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an  extended  Kalman  filter-linked  together  by  a  set  of 
statistical  switching  rules.  The  underlying  idea  of  the 
algorithm  is  to  use  the  batch  filter  to  provide  good  initial 
estimates  and,  when  needed,  reinitialization,  and  to  use  the 
sequential  filter  to  do  the  bulk  of  the  tracking.  The  switching 
rules  are  intended,  in  the  batch  filter,  to  determine  as  quickly 
as  possible  when  the  filter  has  converged  to  a  good  solution  and 
switch  to  the  sequential,  while  in  the  sequential  portion  they 
are  intended  to  determine  when  the  sequential  has  lost  track 
and  switch  to  the  batch  for  reinitialization.  Studies  performed 
at  Tracor  have  shown  that,  in  most  cases,  if  a  sequential  filter 
of  the  type  used  in  the  hybrid  is  initialized  with  reasonably 
good  starting  values,  it  will  be  able  to  maintain  even  the  most 
difficult  track.  Figure  2.2.1  contains  a  block  diagram  showing 
the  relationship  between  the  batch  filter,  the  sequential 
filter,  and  the  switching  rules. 

2.2.2  Aspects  of  Sequential  Filter 

Theory.  The  sequential  filter  used  in  the 
hybrid  algorithm  is  a  variant  of  the  extended  Kalman  filter. 

The  target's  time  rate  of  change  of  the  acceleration  vector  is 
modeled  as  a  Gaussian  white  noise  process.  This  equation  has 
the  form 


da  „ 

at  " 


n 


where  n  is  Gaussian  white  noise  with  zero  mean  and  covariance 
qxx  and  qyy  The  solution,  in  the  mean,  to  this  stochastic 
differential  equation  is  given  by  the  following  linear  equations 
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FIGURE  2.2.1 


HYBRID  ALGORITHM  MACRO  STRUCTURE 
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r (t)  =  r(t8)  +  v(t0)At  +  %a(t0 ) At2 
v(t)  =  v(t0)  +  a (t  0 )A t 
a(t)  =  a(t0) 

where  r  and  v  are  respectively  the  position  and  velocity 
vectors.  The  stochastic  derivative  of  any  tracking  frequency 
is, 


where  e  is  Gaussian  white  noise  with  zero  mean  and  covariance 
qf^.  The  solution,  in  the  mean,  to  this  equation  is 


f  -  f 


Thus,  at  any  time,  t^,  the  state  vector  for  the  model  is: 


X(tk) 


*x<tk) 

ry(tk) 

vx<ck> 

vy(tk) 

ay(tk) 


1 
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where  rx,  fy  are  the  x  and  y  position  coordinates;  vx,  Vy  are 
the  x  and  y  velocities;  ax,  ay  are  the  x  and  y  accelerations; 
and  f01,...,f0  are  the  different  center  frequencies  being 
monitored.  The  state  transition  matrix  for  this  system  is 


0 

At 

0 

At2  /  2 

0 

0 

...  0 

1 

0 

At 

0 

At2  /  2 

0 

.  .  .  0 

0 

1 

0 

At 

0 

0 

...  0 

0 

0 

1 

0 

At 

0 

...  0 

0 

0 

0 

1 

0 

0 

...  0 

0 

0 

0 

0 

1 

0 

.  .  .  0 

0 

0 

0 

0 

0 

1 

...  0 

0 

0 

0 

6 

0 

0 

...  1 

where  R?X6  is  a  matrix  containing  transition  parameters  for 
position,  velocity,  and  acceleration  only,  0Sxn  and  0nx  are 
matrices  of  zeroes,  and  Inxn  is  the  identity  matrix.  Note  that 
this  corresponds  to  the  constant  acceleration  motion  model  used 
in  the  MLP.  However,  the  time  periods  over  which  the  transition 
matrix  is  applied  are  far  shorter  than  the  time  periods  over 
which  the  MLP  maneuver  models  are  applied.  This  is  one  of  the 
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major  advantages  of  the  sequential  filter  over  the  MLP.  There 
is  no  need  to  test  for  appropriateness  of  model. 

The  measurement  models  are  the  same  as  those 
used  in  the  MLP.  Again,  the  primary  models  are  Doppler-shifted 
frequency  and  bearing,  with  the  capability  to  use  range,  Doppler 
ratio,  Doppler  difference,  and  time  difference  of  arrival. 

In  order  to  obtain  a  practical  propagation  and 
update  algorithm,  the  extended  Kalman  filter  uses  a  linearized 
version  of  the  measurement  model,  i.e.,  a  first  order  Taylor 
series  expanded  about  the  most  recent  state  estimate.  For 
example,  suppose  a  measurement  of  Doppler-shifted  frequency, 
fr^ ,  was  received,  then: 

^rj  =  ^rj  (tk-i  ’  ^ck-i^ 

3  f  .  ~ 

+  -5“  (ck.i  ■ 
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where , 


tk  *  time  of  current  estimate 


frj (t^jX(t^))  =  observed  Doppler-shifted  value 

of  frequency  r j 


3f. 


3(0 


tk-1  =  time  of  last  estimate 

—  (tk  ,X(tjc_i )}  =  partial  derivative  of  received 

Doppler-shifted  frequency  frj  with 
respect  to  (•)  evaluated  at 

Ck-i  ’^(tk-i ) 


6^  =  Knonecker  delta  = 


1  if  i“j 
0  if  i^j 


Bearing  measurements  can  also  be  expanded  in  a  first  order 
Taylor  series  containing  only  position  terms.  Letting 


dk  -  frj(tk,X(tk))  -  frj(tk.,,X (tk.,)) 


we  have,  in  vector  form: 

A 

dk  =  Hk(x<tk)  -  X(tk_1)) 

where , 
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Hk  =  (  ax'*  (Ck-1  ,X^ck-i^ . 


3fr 

3a, 


i(tk-.’5<tk-.>)' 0 . 


3^ri  /  ■+■  \ 

3 f0^  rk-i  ’^k-iV’  •  •  •  ,0/ 


X(Ck)-X(tk_1)  =  ^x(tk)-x(tk_1)  ,  .  .  .  ,a(tk)-a(tk_1)  ,0, 

^oj  (ck) ~^°j  ^ck-i ) » • • • ^ ) 

The  error  covariance  matrix,  P(tk) ,  for  the 
parameters  is  given  by: 

P(tk)  -  *(tk,tk.1)P(tk.1)*T(tkltk.I)  +  r(tk,tk_3) . 


The  matrix,  r(tk,tk_1),  contains  the  effects  on  the  covariance 
of  process  noise  and  is  found  by: 


P(ck»ck_i)  = 


<l,(t,T)Q(T)<j)T(t,T)dT. 


This  equation  can  be  solved  analytically  because  Q  is  assumed 
to  be  a  constant  matrix  of  the  form: 
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Figure  2.2.2  contains  a  logical  flow  diagram  of 
the  sequential  processor. 

Implementation .  There  are  several  schemes  which 
may  be  used  to  implement  the  sequential  filter,  each  with  its 
own  set  of  advantages.  The  method  of  implementation  chosen  for 
the  hybrid  algorithm  is  the  UDUT  square  root  form  of  the  filter. 
For  this  particular  problem  it  possesses  several  advantages: 

(1)  Helps  prevent  conditioning  problems  by 
insuring  that  the  covariance  matrix,  P, 
will  be  positive  definite  and  symmetric 
at  all  times. 

(2)  Helps  insulate  algorithm  from  effects  of 
changes  in  computer  word  size,  making  it 
more  suitable  for  operation  on  minicomputers 
than  other  forms  of  the  algorithm. 
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xi<->  - 


Estimated  state  vector  at  time  t^  before  measurement 


xt(+) 


Estimated  state  vector  at  time  t^  after  measurement 


ci 


=  Time  t- 


Pi(-) 


Estimated  covariance  matrix  at  time  t^  before  measurement 


f! 

•  i 


P^(+)  =  Estimated  covariance  matrix  at  time  t^_  after  measurement 


Zi 


Di 


=  Measurement  (Doppler  shifted  frequency,  bearing,  etc.) 


at  time  ^ 


Measurement  residual  at  time  t^ 


FIG.  2.2.2 


LOGIC  DIAGRAM  -  KALMAN  FILTER 


28 


v 


Tracor  Applied  Sciences 


(3)  Minimizes  storage  requirements. 

(4)  Is  compatible  with  Givens'  rotation  form 
of  MLP . 

A  summary  of  the  standard  Kalman  Filter  conditions 

and  equations  is  given  in  Figure  2.2.3.  The  following  will 

T 

describe  how  these  equations  are  implemented  under  the  UDU 
decomposition . 


defined  as 


Recall  that  in  Section  2.1.2  the  covariance  is 


where  U  is  upper  triangular  with  ones  on  the  diagonal  and  D  is 
diagonal.  This  gives  rise  to  the  following  definitions  for  the 
sequential : 


Pk<-)  - 

pk(+)  =  u1;1(+)d*1(+)u-tc+) 

Note  the  following  from  Figure  2.2.3: 

(1)  H^(-)  is  a  lxn  vector. 

(2)  R^  is  a  scalar. 

(3)  Q  is  constant  for  all  t. 

Substituting  into  the  expression  for  the  Kalman  gain  we  get: 
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System  Model: 

Xk  =  Ktk.tk_i)Xk_i  +  w(t)  w(t)  ~  N(0,Q) 

Measurement  Model: 

zk  ■  \(XM  +  \  \  -  N<°'Rk> 

State  Estimate  Propagation: 

5k(-)  -  ®Cck.ck-l>xk-l<+> 

Error  Covariance  Propagation: 

f  tk 

pk<-)  -  *(tk,Ck.1>pk.1(+)»Tftk,tk.1)  +  ;tk  i»(c,T) q(t)»t< 

Update  estimates  at  time  =  tk  based  on  measurement  Zk. 

State:  \(+)  =  Xk<">  +  Kk[Zk  ~  hk(xkC-))  3 

Covariance:  Pk(+)  =  pk<“)  “  KkHk(-)pk(-) 

Gain  Matrix:  Kk  =  Pk(-)Hj(-) tHk(-)Pk(-)Hjj:(-)  +  Rk] _1 

where  hk  =  measurement  model  used  at  time  k  (Doppler  shift, 
bearing,  etc.) 

x<tk>  =  xk(-) 

FIG.  2.2.3  -  KALMAN  FILTER  EQUATIONS 


_  5hk<x(tk>) 

Hk(_)  “  3X(tk) 


,t)dt 
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I 


m 


-  ur1 (-)Dr1 (-)u‘t(-)h't(-> 


Hk(’>uk1<->Dk1(-)ukT(-)HkT(->  + 


Let  vj(-)  =  H,  (-)Ur1  (-)  , 


-l 


\T' 


Kk  =  U-(-)D-'(-)v-T(-) 


vJ(-)dCx(-)v.  (-)  + 


Since  v£(-)  is  l*n,  D^(-)  is  n*n ,  and  is  a  scalar,  we  have: 

a  =  V‘T(-)D£X(-)V£l(-)  +  (a  =  a  scalar). 


Thus , 


Kk  - 


uk-'(-)^‘(-)vk(-) 


Substituting  this  expression  into  the  state  estimate  update 
gives : 


Xk(+)  -  Xk(-)  + 


Substituting  into  the  covariance  update  gives 


Pk(+)  -  Pkc->  -  Vk(-)Pk(-) 
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uk(+)Dk1(+)ukT(')  =  uk1(")Dk1(")ukT(_) 


uk'  <-)\(-)Hk(-)«k'  <-)»;’  <-)“kT<-> 


Dk1(-)\1(')Vk(-)Dk(') 

D  (-)  -  -£ - £ - £ - £ - 

k  v  a 


\V)  (2) 


For  the  error  covariance  propagation  the 
decomposition  gives: 

PkW  '  *(tk,tk-l)Pk-i(+,‘fT(tk’tk-i)  +  r<tk,tk-l) 


♦<tk>v.>’,il.(+>C.<+>C«+Ml(tk>tk-.> +  r(tk-tk-.> 


Let 


sk-i(+>  -  »<tk'tk-.)0kl!(+)-  then 


pk<->  •°k1<-)Dk  ■  sk-,(+)Dk-l(+)sk-.<+)  +r<tk-tk-.)  <3’ 


Given  the  standard  Kalman  Filter  equations 
expressed  in  square  root  format,  the  procedure  is  to  solve 
equation  (3)  for  Uk(-)  and  Dk(-).  This  is  an  easy  and  straight¬ 
forward  procedure  to  carry  out  since : 

(1)  SiL  =  1  because  4’ii(tk,tk_1 )  - 

Pku(->  =  Pkii(+)  =  1;  where  is 

an  element  in  U”1  . 
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(2)  The  linear  equations  generated  by 
equation  (3)  are  all  solvable  by  back 
substitution. 

(3)  For  j <i ,  ^kii<'>  =  for  i=J .  ^k^C")  = 

T 

(-)  =  1,  allowing  all  the  d^^C-)  terms 
to  be  found.  For  j>i,  can  be  found 

since  all  the  dk^(-)  terms  required  are 
known. 

Once  the  U^1 (-)  and  D^1 (-)  matrices  have  been  determined,  they 
can  be  substituted  in  equation  (1)  to  give  the  updated  estimate 
and  substituted  into  equation  (2)  to  give  a  set  of  equations  to 
be  solved  for  Uj^1  (+)  and  D^.1  (+)  that  use  a  method  identical  to 
the  one  employed  to  get  U^1 (-)  and  D^1 (-) . 


are : 


Other  considerations  in  implementing  the  filter 


(1)  The  matrix  <J|(t^,t^_i)  is  never  stored. 
Whenever  a  multiplication  is  required 
the  special  nature  of  3>  is  exploited. 

(2)  The  V  (t^'t^.j )  integration  is  implemented 
in  a  closed  form  solution  requiring  no 
numerical  integration  techniques. 

(3)  Storage  requirements  are  kept  to  a 
minimum.  D^1  is  stored  as  vector  and 
U^1  is  stored  in  row  symmetric  form 
without  the  diagonal  elements. 
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2.2.3  Aspects  of  the  Batch  Filter 

Theory.  The  batch  portion  of  the  hybrid 
algorithm  is  based  on  the  assumption  that,  for  short  periods  of 
time,  the  motion  of  a  submarine  can  be  approximated  by  the 
constant  acceleration  model  of  Section  2.2.2  (this  model  is 
quadratic  in  t,  allowing  it  to  fit  tracks  which  are  moderately 
curved) .  Once  it  has  been  determined  that  the  batch  algorithm 
has  converged  to  a  good  solution,  a  switch  is  made  to  the 
sequential  filter.  Thus,  in  the  hybrid  algorithm,  the  batch 
filter's  job  is  to  converge,  as  quickly  as  possible,  to  a  state 
vector  which  will  allow  the  sequential  filter  to  be  properly 
initialized. 


The  measurement  models  used  in  the  batch  filter 
are  the  same  as  those  used  in  the  MLP  and  in  the  sequential 
filter,  mainly  bearing  and  Doppler-shifted  frequency.  As  in 
the  MLP,  the  batch  filter  uses  data  measurements  to  provide 
state  vector  estimates.  The  measurement  models  are  approximated 
by  first  order  Taylor  series  expansions.  The  design  (X)  matrix 
is  made  up  of  the  corresponding  partial  derivatives,  and  the 
variation  in  the  state  vector  is  the  quantity  estimated  (see 
Section  2.1.1).  As  in  the  MLP,  a  major  share  of  the 
computational  effort  involved  in  the  filter  is  in  inverting 
the  X  a-1X  matrix,  and,  as  in  the  MLP,  this  is  done  using  a 
sequence  of  Givens'  rotations  (Section  2.1.1).  In  addition  to 
the  benefits  described  in  Section  2.1.1,  there  is  one  further 
and  very  substantial  benefit  to  using  this  method  in  the  hybrid 
algorithm.  Specifically,  the  U_1  matrix  of  the  sequential  is 
the  inverse  of  the  U  matrix  in  the  batch  filter,  and  the  D-1 
matrix  of  the  sequential  is  the  inverse  of  the  batch  D  matrix 
(Sections  2.1.1  and  2.2.1),  i.e.. 
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U_1  =  TT-1 

U(seq)  u (batch) 
D(seq)  ~  D(batch) 


are  unit  upper  triangular  and  that  D(seq) 


Note  that  U  and  U 
and  D (batch)  are  diagonal,  thus  their  inverses  are  easily  found 
when  switching  from  the  batch  to  the  sequential . 


Since  the  batch  filter  contains  only  one  motion 
model,  there  are  no  tests  needed  to  select  between  various 
models  and  no  tests  needed  to  select  a  data  interval  over  which 
a  given  model  is  appropriate.  Since  the  motion  model  can  be 
expressed  in  terms  of  the  transition  matrix  used  in  the 
sequential  filter,  all  forward  and  backward  mappings  can  be 
performed  with  the  same  code  used  in  the  sequential  filter. 
These  two  features  reduce  substantially  the  computer  overhead 
imposed  by  the  batch  filter  as  compared  with  the  overhead 
imposed  by  the  MLP .  Figure  2.2.4  contains  a  flowchart  of  the 
logic  flow  in  the  batch  filter. 

Implementation .  As  discussed  above,  several  of 
the  computational  aspects  of  the  batch  filter  are  identical  to 
those  in  either  the  sequential  filter  or  the  MLP.  Briefly, 
these  are: 


(1)  Solving  the  approximating  linear  least 
squares  problem  using  Givens'  rotations 
(MLP,  Section  2.1.1). 
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FIG.  2.2.4  -  FLOWCHART  -  HYBRID  BATCH  INITIALIZER 
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(2)  Projecting  the  track  using  a  constant 
acceleration  model  (sequential,  Section 
2.2.1)  . 

(3)  Computing  partial  derivatives  for  a  given 
measurement  model  (MLP,  Section  2.1.1). 

For  a  given  number  of  data  points,  the  batch 
filter  generates  a  sequence  of  initial  states,  stopping  when 
the  sequence  converges  to  a  particular  state  vector.  The 
sequence  is  said  to  have  converged  when: 


XQH  ~ 

X«H 


£  E  for  j  =  1,2,  ...  ,n 


where , 


=  the  j—  component  of  the  state  vector 
after  the  i—  iteration. 

If  the  sequence  hasn't  converged  after  30  iterations  (i>j),  a 
new  point  is  read  and  the  process  is  begun  again. 

One  significant  computational  difference  between 
the  MLP  and  the  batch  filter  is  in  minimizing  the  residual  sum 
of  squares  once  the  variation  in  X  has  been  found,  using  weighted 
least  squares.  The  MLP  uses  the  Marquardt  (Reference  5) 
algorithm  to  speed  convergence ,  while  the  batch  filter  uses  a 
one -dimen s ional ,  quadratic  search  procedure.  The  algorithm  is 
extremely  fast  and  has  several  important  properties: 
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(1)  It  searches  in  one  direction  only,  avoiding 
the  complications  of  multi-dimensional 
searches . 

(2)  For  a  quadratic  function,  the  method 
converges  quadratically  to  the  minimum. 
Since  even  non-quadratic  functions  behave 
approximately  quadratic  in  the  region  of 
a  minimum,  this  assures  rapid  convergence 
in  the  final  states  of  computation. 

(3)  It  requires  only  information  at  the 
present  stage  and  the  one  immediately 
previous  to  the  present,  thus  reducing 
substantially  the  amount  of  core  and 
computation  required. 

(4)  The  method  requires  no  derivatives. 

2.2.4  Switching  Rules 

2 . 2 . 4 . 1  Batch  to  Sequential.  In  the  hybrid 

algorithm  the  basic  assumptions  are  that  a  simplified  MLP 
algorithm  can  be  used  to  fit  the  data,  the  fit  can  be  detected 
early,  and  a  switch  can  be  made  to  a  sequential  filter  using 
the  initial  conditions  supplied  by  the  batch  filter.  Thus, 
the  batch-to-sequential  switching  must  try  to  satisfy  two 
constraints : 


(1)  The  transition  from  the  batch  filter  to 
the  sequential  filter  should  be  made  as 
quickly  as  possible. 
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(2)  The  switch  from  batch  to  sequential 

should  not  be  made  until  a  good  approxi¬ 
mation  to  the  initial  state  vector  has 
been  found. 

To  resolve  these  problems  the  hybrid  algorithm 
uses  a  nonlinear  test  of  significance  developed  by  Gallant 
(Reference  6) .  The  Gallant  method  is  used  to  test  the  hypothesis 

H0  :  X^  =  Xj  versus  Hj  :  X^_  ?  Xj 

where  X^  is  the  estimated  state  vector  from  the  first  i  data 
points  and  Xj  is  the  estimated  state  vector  from  the  first  j 
data  points.  Currently,  i  =  j  +  2*NBY,  where  NBY  is  the  number 
of  buoys  in  the  scenario.  Gallant's  test  statistic  for  this 
hypothesis  is: 

j'k  •  2 

T  _  k^l _ 

>  j  )  j  I  -+  \  , 

I,»k  •  (^k-VXi)) 1 
k=l 

where 

Yv:  =  k—  observed  measurement  value 

~  -*■  t;U 

Y^(Xp)  =  predicted  k —  measurement  value  based 
on  estimated  state  vector  Xp(p=i,j) 

t"  H 

=  weight  associated  with  k—  measurement 
The  critical  value  for  the  test  is: 
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where 


C(i)  -  1.  +  l»Fa(P.i--P) 

1-p 


p  =  number  of  elements  in  the  state 


vector 


i  =  number  of  points  used  to  estimate 
state  vector  X,- 


Fa(p,i-p)  =  a  percentage  point  of  F-distribution 


Lastly,  the  decision  rule  is: 


(1)  If  T  is  greater  than  C,  reject  H0 

(2)  If  T  is  less  than  C,  accept  H0 . 


i  I 


Basically,  the  idea  is  to  compare  the  weighted 
residual  sum  of  squares  for  estimates.  If  the  ratio  is  too 
large,  we  reject  the  idea  that  the  estimates  are  equal,  otherwise, 
we  accept  it.  Gallant's  critical  point  gives  us  a  convenient  and 
consistent  measure  of  "too  large". 


If  i  =  j  +  2*n  and  k  =  i  +  2*n,  where  n  is  the 
number  of  buoys,  then  a  switch  is  made  from  the  batch  to  the 
sequential  when 


T(i . j )  -  C(i)  and  T(k,i)  1  C(k)  • 


!  I 
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Figure  2. 2. 4.1  is  a  flowchart  of  the  batch-to- 
sequential  switching  logic.  There  are  two  primary  reasons  for 
requiring  the  data  to  satisfy  the  switching  rules  two  times  in 
succession: 

(1)  By  chance ,  a  certain  percentage  of  the 
time  the  data  will  give  a  false  reading 
due  to  the  inherent  random  nature  of  the 
data . 

(2)  During  the  earliest  stages  of  estimation 
estimates  of  the  state  vector  may  be  quite 
poor  and  may  not  change  much  over  2*n 
points.  Thus,  based  on  the  data,  the  test 
may  conclude  that  the  batch  filter  has 
converged  when,  in  fact,  it  hasn't. 

Usually  another  run  through  2-n  points 

is  enough  to  determine  whether  or  not 
the  process  has  converged. 

The  basic  idea  behind  this  approach  is  to  use  the  data  itself 
to  determine  when  the  batch  filter  has  converged.  At  present 
an  alpha  level  of  .005  is  being  used,  however,  this  could 
become  an  operator  input  if  desired. 

2. 2. 4. 2  Sequential  to  Batch.  The  sequential  filter 

used  in  the  hybrid  algorithm  is  capable  of  maintaining  track 
under  a  wide  range  of  conditions.  However,  studies  at  Tracor 
have  shown  that  there  are  times  when  it  will  lose  track,  for 
example,  if  there  is  a  long  period  of  time  when  data  is  received 
from  only  one  sensor  due  to  low  signal-to-noise  ratios.  Thus, 
the  job  of  the  sequential -to-batch  switching  rules  is  to  sense 
when  the  sequential  filter  has  lost  track  and  cause  a  transfer 
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FIG.  2. 2.4. I  -  FLOWCHART  -  BATCH  SWITCHING  LOGIC 
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back  to  the  batch  filter  for  reinitialization.  Note,  an  implied 
constraint  is  that  the  procedure  must  not  be  so  sensitive  as  to 
interpret  normal  maneuvers  as  losing  track. 

The  hybrid  contains  two  tests  to  determine  when 
the  sequential  filter  has  lost  track.  The  first  is  a  gross 
test,  and  failure  to  pass  it  will  cause  an  immediate  switch 
to  the  batch  filter.  Essentially,  it  consists  of  a  comparison 
of  the  residuals  obtained  before  and  after  the  filtering.  If 
the  residual  after  correction  is  larger  than  the  residual  before 
correction,  the  filter  is  assumed  to  be  diverging  and  a  switch 
is  made  to  the  batch.  This  is  a  fairly  insensitive  test; 
however,  it  is  fast,  easy  to  do,  and  provides  "worst  case" 
protection. 


The  second  test  has  proved  to  be  quite 
successful  in  determining  loss  of  track.  It  is  based  on  two 
aspects  of  the  Kalman  Filter: 

A 

(1)  If  state  vector  estimate  X^_1  is  unbiased, 
then  state  vector  estimate  is  unbiased 
(Reference  7),  i.e.,  E[^]  =  E [ ]  . 

(2)  The  errors  associated  with  the  measurement 
model  are  assumed  to  be  normally  distributed 
with  zero  mean. 

till 

Define  the  i—  measurement  residual  to  be 

Yi  -  Yi«i) 
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where 


th 

=  i—  observed  measurement  value 

Y^(X^)  =  i—  predicted  measurement  value  based 
on  estimated  state  vector  X^ 

The  residuals  are  the  differences  between  what  was  observed  at 
the  i—  data  point  and  what  the  filter  predicted^at  that  point, 
and  may  be  thought  of  as  the  observed  errors  if  X^  is  near  X^. 
Under  aspect  (2)  above,  the  residuals  should  be  distributed  as 
normal  random  variables  with  zero  mean  if,  again,  Xfc.  is  near  X^. 
Specifically,  if  the  residuals  are  plotted  over  time  and  exhibit 
no  trends,  such  as  ramping  or  dramatic  mean  shifts,  then  X^  is 
probably  near  Xk;  at  least  the  data  does  not  violate  this 
assumption. 


Thus,  the  problem  of  sequential-to-batch 
switching  may  be  viewed  as  a  problem  in  determining  whether  or 
not  a  trend  exists  in  the  measurement  residuals.  To  do  this 
the  hybrid  fits  a  regression  line  of  the  form, 

y=b0+b1t  (t  =  time) 

to  the  residuals  from  each  sensor  for  each  data  type.  The 
lines  are  tested  for  significance  using  an  F-test  (currently 
a  equals  .005)  and,  if  any  one  of  them  tests  as  significant  for 
four  consecutive  data  points  a  switch  is  made  to  the  batch  filter. 

Note  that  a  separate  line  is  computed  for  each 
data  type  from  each  sensor.  This  is  done  to  prevent  information 
from  one  or  two  sensors,  which  may  be  sensitive  to  the  filter 
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losing  track,  from  being  buried  or  obscured  by  the  information 
from  several  other  sensors  which  are  insensitive  to  the  filter 
losing  track  due  to  scenario  geometry,  poor  signal-to-noise 
ratios,  etc.  There  are  several  reasons  for  requiring  the  line 
to  test  significantly  for  more  than  one  data  point: 

(1)  Early  in  the  data  stream  it  is  very 
possible  that  the  residuals  will  cause 
a  significant  fit  to  occur  when,  in 
fact,  track  has  not  been  lost. 

(2)  The  sequential  filter  will  typically 
wander  off,  at  different  times,  for 
one  or  two  points.  This  allows  it 
time  to  regain  track  on  its  own. 

(3)  One  or  two  outliers  at  any  point  in  the 
data  stream  could  cause  an  indication 
of  significance  when,  in  fact,  there  is 
none . 

(4)  It  prevents  sharp  or  rapid  maneuvers 
from  being  interpreted  as  lost  track. 

To  allow  the  hybrid  to  be  altered  to  fit 
various  environmental  conditions,  it  would  be  quite  easy  to 
allow  the  F-test  significance  level  and  the  number  of  required 
significant  points  to  be  operator  inputs.  Then,  for  example, 
if  the  environment  was  noisy  with  low  data  rate  the  alpha  level 
could  be  reduced  and  the  number  of  significance  points  could 
also  be  reduced,  causing  the  sequential  to  be  more  sensitive 
to  any  loss  of  track. 
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2.3  Iterated  Sequential  Algorithm 

2.3.1  Aspects  of  Iterated  Sequential  Algorithm 

The  iterated  sequential  algorithm,  in  form, 
is  very  similar  to  the  hybrid  algorithm.  That  is,  both 
algorithms  are  comprised  of  three  primary  elements--an 
initializer,  a  sequential  Kalman  Filter,  and  a  set  of  switching 
rules  for  moving  from  the  initializer  to  the  tracker  and  back 
again.  The  only  difference  in  the  two  procedures  is  in  the 
initializer.  The  hybrid  uses  a  batch  algorithm  to  initialize, 
while  the  iterated  sequential  uses  the  same  filter  to  initialize 
that  it  does  to  track,  relying  on  iteration  to  produce 
convergence.  There  were  several  reasons  for  implementing  this 
algorithm: 


(1)  There  is  a  reduction  of  core  requirements. 
Instruction  requirements  are  reduced 
because  the  same  subroutines  used  in  the 
straight  sequential  filter  are  used  in 
the  initializer.  Data  requirements  are 
reduced  because  fewer  save  arrays  are 
needed  to  make  the  transition  between 
initializer  to  tracker. 

(2)  There  is  faster  execution  time.  The 
iterated  sequential  starter  uses  fewer 
operations  to  achieve  an  estimate  than 
does  the  batch,  however,  iterating  about 
an  answer  could  negate  this  to  some 
degree . 
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(3)  One  possible  weakness  of  the  hybrid  is 
attempted  initialization  during  a  track 
too  difficult  to  follow  with  the  constant 
acceleration  model.  By  its  sequential 
structure  the  iterated  starter  should  be 
able  to  avoid  this  problem. 

(4)  This  algorithm  has  simpler  program 
organization.  In  the  hybrid,  provisions 
must  be  made  to  accommodate  two  distinctly 
different  algorithms  and  transfers  between 
them.  In  the  iterated  sequential,  the 
same  operations  are  performed  in  both 
portions  of  the  algorithm,  simplifying 
considerably  storage  management  and 
control  transfer  between  the  initializer 
and  tracker. 

Figure  2.2.1  can  equally  well  represent  the 
iterated  sequential  algorithm  as  the  batch;  simply  replace 
"BATCH"  with  "ITERATED  SEQUENTIAL".  The  extended  Kalman  Filter 
used  for  tracking  is  identical  to  that  used  in  the  hybrid  (see 
Section  2.2.2)  and  the  switching  rules  (Section  2.2.4)  are  also 
the  same . 


Figure  2.3.1  contains  a  logical  flowchart  of  the 
iterated  sequential  starter.  The  basic  idea  is  to  use 
accumulated  data  points  to  project  initial  state  vector  and 
covariance  estimates  forward,  using  the  sequential  filter, 
and  then  test  for  convergence.  Convergence  occurs  when 


x(j)i  -  xq)i-i 
XCj)t 


E  for  j  =  1 ,2 ,3  , .  .  .  ,n 
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FLOWCHART  -  ITERATED  SEQUENTIAL  INITIALIZER 
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where , 

cli 

X(j)i  =  the  j  —  component  of  the  state  vector 
after  the  i—  iteration 

If  the  convergence  does  not  occur  on  the  i—  iteration,  the 
state  vector  is  projected  back  to  time,  t0 ,  using  the  constant 
acceleration  model  (Sections  2.1  and  2.2.2)  as  implemented  in 
the  sequential  filter.  When  convergence  does  occur,  a  switch 
test  is  made.  If  the  decision  is  made  to  switch,  the  straight 
sequential  filter  is  entered;  if  not,  the  state  is  projected 
back  (as  above),  a  new  data  point  is  read,  and  the  cycle  begins 
anew. 
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3.0  TEST  RESULTS  (COMPARISONS) 

3 . 1  Introduction 


This  section  presents  the  results  of  several 
tests  designed  to  measure  the  relative  tracking  abilities  of 
each  of  the  three  algorithms  examined  in  this  report.  The 
tests  consisted  of  four  different  scenarios,  each  replicated 
twenty-five  times  as  a  series  of  Monte  Carlo  runs.  The 
algorithms  then  processed  each  scenario  and  the  results  were 
analyzed  according  to  four  criteria.  These  requirements  are: 

(1)  Average  distance  error,  <T,  throughout 
the  track,  defined  as, 

5  =  (t£  S_(t:)dt 

K.  tf -ti 

where , 

5(t)  =  rE(t)  -  rT(t) 

rjr(t)  =  estimated  position  vector  at  time,  t 
rT(t)  *  true  position  vector  at  time,  t 

This  is  a  measure  of  the  average  closeness 
of  the  true  and  estimated  trajectories  for 
the  time  interval  (tj[,tf).  The  two 
trajectories  are  considered  close  when 

T  <  6* 
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For  some  specified  level  of  6*  (currently 
500  meters)  . 

(2)  Percent  holding  time,  ”,  throughout  the 
track,  defined  as, 


.1  (Tfj-Ti .) 

T  =  — - -  •  100% 

zf  ~zi 

where,  for  a  given  time  interval  (t^,tf), 
(T ij  ,Tfj)  C  (ti,tf) 


(Tij  .ffj)  n  (Tik,tfk)  =  <f>  for  j  +  k 


6j  =  /fj  4 '-rdt  1  5*  for  j  =  1-2 . n- 

J  7ij  fJ 


(3)  Predictive  ability.  Assuming  the  target 
stays  on  the  same  trajectory,  a  given 
algorithm  will  be  said  to  have  good 
predictive  ability  if 

<5(tf+At)  <_  5*. 

(4)  Average  CPU  time  used  to  process  one 
scenario.  This  gives  a  measure  of  the 
relative  speeds  of  the  three  algorithms. 


Note,  again,  that  these  measures  are  computed  over  twenty-five 
Monte  Carlo  samples  and  are  intended  as  a  measure  of  the  average 
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or  expected  tracking  abilities  of  the  algorithms. 

3 . 2  Scenario  Descriptions 

In  testing  the  MLP,  hybrid,  and  iterated 
sequential  algorithms  the  intent  was  to  produce  several 
scenarios,  covering  a  range  of  tracking  difficulties,  to  test 
the  capabilities  of  the  three  algorithms.  The  intent  was  not, 
at  this  time,  to  evaluate  each  algorithm  at  the  extreme  edge  of 
its  performance  range . 

The  basic  layout  for  each  of  the  four  scenarios 
generated  was  the  same.  It  consisted  of  three  DIFAR  buoys 
arranged  in  the  shape  of  an  equilateral  triangle  with  2,000 
meter  sides;  no  buoy  drift  was  assumed.  All  maneuvers  were 
assumed  to  take  place  within,  or  near,  the  buoy  field  and  a  sea 
state  of  two  was  assumed.  Figure  3.2e  summarizes  the  environ¬ 
mental  conditions  and  inputs  common  to  all  four  scenarios. 

The  first  scenario  run  (Figure  3.2a)  was  a 
simple  straight  line  course  through  the  center  of  the  buoy 
field.  The  intent  of  this  scenario  was  to  determine  whether  or 
not  the  three  algorithms  could  track  straight  line  motion  with 
high  efficiency. 

The  second  scenario  considered  (Figure  3.2b)  was 
a  180  degree  turn  in  the  middle  of  the  buoy  field.  This  was 
felt  to  be  a  scenario  of  moderate  difficulty  and  was  intended 
to  compare  the  ability  of  each  algorithm  to  respond  to  an 
extremely  sharp  turn.  Note  that  the  maneuver  time  for  this 
turn  is  about  300  seconds.  This  not  only  includes  the  time 
required  to  make  the  turn  itself,  but  also  the  time  needed  to 
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regain  the  required  operating  velocity  (in  this  case  7  meters/ 
second) . 


Scenario  3  (Figure  3.2c)  was  considered  to  be 
of  high  difficulty  and  consisted  of  a  360  degree  turn  performed 
inside  the  buoy  field.  This  was  intended  to  test  the  ability 
of  each  algorithm  to  follow  a  severe  maneuver  for  a  prolonged 
length  of  time.  Again,  the  maneuver  time  indicated  in  the 
figure  included  not  only  the  time  required  to  make  the  circle, 
but  also  the  time  required  to  regain  an  operating  velocity  of 
7  meters/ second. 

The  final  scenario  (Figure  3. 2d)  was  also 
considered  to  be  of  severe  difficulty  and  was  intended  to 
measure  the  ability  of  each  algorithm  to  follow  a  sequence  of 
fairly  rapid  turns  and  also  to  compare  the  ability  of  each 
algorithm  to  follow  maneuvers  outside,  but  near,  the  buoy  field. 

3.3  Generation  of  Data  for  DIFAR  Simulation 


In  order  to  carry  out  simulations  of  the  three 
target  tracking  algorithms  examined  in  this  report,  an  algorithm 
for  generating  realistic  data,  as  it  would  appear  from  a  general 
processor  for  DIFAR  buoys,  was  developed.  It  was  desired  that 
the  generated  data  {  (f^.B-^.t^)}  have  a  statistical  distribution 
closely  approximating  real  data. 

3.3.1  General  Approach 

The  general  processor  modeled  (see  Figure  3. 3. 2.1) 
consists  of  a  filter  band  with  bins  of  width  Af ,  followed  by  a 
square  law  detector  [(  )2  and  a  1/Af  averager]  and  a  finite 
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FIG.  3.2b  -  SCENARIO  2.  DESCRIPTION 


Duration:  765  seconds 
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FIG.  3.2c  -  SCENARIO  3  DESCRIPTION 
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FIG.  3. 2d  -  SCENARIO  4  DESCRIPTION 
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COMMON  SCENARIO  INPUT  AND  ENVIRONMENTAL  VARIABLES 


FIGURE  3. 3. 2. I  -  GENERAL  PROCESSOR  FOR  DIFAR  BUOY 
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time  perfect  integrator  (ALI) .  A  MAX-OR  picks  the  strongest 
signal  in  frequency  from  the  omni  signal  channel  and  computes 
bearing  from  the  X  and  Y  channels.  The  following  describes 
the  statistics  of  the  outputs  of  various  components  in  such  a 
processor.  The  distribution  of  the  envelope  of  a  Gaussian 
process:  (1)  at  the  output  of  a  narrowband  filter  is  Rayleigh 

distributed,  (2)  at  the  output  of  a  narrowband  filter  and 
square  law  detector  is  exponential,  and  (3)  at  the  output  of 
a  narrowband  filter,  square  law  detector,  post  detection 
integrator  is  CHI -Squared  distributed.  One  could  randomly 
generate  errors  from  these  distributions  to  add  to  the  correct 
data;  however,  with  signal  present,  the  X  and  Y  signals  are 
not  independent . 

3.3.2  Analysis 

3.3. 2.1  Generation  of  Frequency  Data.  It  is 

assumed  that  the  input  to  the  processor  may  be  written  as  the 
sum  of  a  narrowband  Gaussian  signal  and  Gaussian  noise, 

W(t)  =  S(t)  +  n (t) 


where 


S(t)  -  N(0,cs2) 
n(t)  -  N(0,an2) . 

At  the  output  of  a  narrowband  filter,  the  signal 
term  may  be  represented  as  the  sum  of  two  orthogonal  signals, 
each  multiplied  by  an  independent  Gaussian  random  variable. 
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where 


S(t)  =  S  costu0t  +  S2sino)0t 

51  -  N(0,as2) 

52  -  N(0 ,asl) 


0)o  =  center  frequency  of  filter 
Similarly  for  the  noise  term, 


where 


n(t)  =  r^cosoiQt  +  n2sino)0t 


n:  ~  N(0,an2) 
n2  ~  NCO.a^2) . 

In  the  processor,  W(t)  is  squared  and  time 
averaged  over  1/Af  seconds,  where  Af  is  the  filter  bandwidth  in 
Hertz . 
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-  ,1/Af 

W(t)2  =  Af  [S  2cos2aj0t  +  S22sin2cuQt  +  n12sin2ca0t 

•'0 

+  n22cos2ajQt  +  2S1nJcos2<u0t  +  2S2n  2cos2o>0t 
+  2SX  S2coso)0tsinw0t  +  2ri  i  n  2cosoj0tsinoj0t 
+  2S:  Ti2cosa)0tsiruu0t  +  2S2n  l  coscu0tsinaj0t  ]  dt . 

It  is  assumed  that  S2 ,  S2>  ,  and  n2  are  slowly  varying 

compared  to  cos2ojQt,  sin20)ot,  and  coscu0tsinw0t .  These  may  then 
be  pulled  outside  the  integral  sign.  The  terms  with  cosaj0tsinco0t 
factors  integrate  to  zero. 


-  ,1/Af 

W(t)2  =  Af [ ( S 1 2  +  n22  +  2SX  n  2 ) J  cos2wQt 


dt 


1/Af 


r  j. 

+  (S22  +  n j 2  +  2S2n2)J  sin2w0t  dt] 


Integration  yields 


W(t)  2  =  Af  [  (Sj 2  +  n  2  +  2SX  nx )  (%t  +  r—  sin2o>0t) 


+  (S22  +  nt2  +  2S2n2)  (%t  -  z~~  sin2u0t)] 


1/Af 

0 


Since  the  process  is  narrowband,  Af  <<  coQ 


W(t)2  =  % ( S2  2  +  n22  +  2Sir'i) 


+  %(S2  2  +  rij 2  +  2S2n  2) 
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Using  the  transformations 


S  *  -  o  gN  3 
S2  =  aSN2 

Hi  *  %N3 
n2  =  a^N,, 


where  N3  ,  N2 ,  N3  ,  N4 


N(0,1),  the  expression  becomes 


W(t)2 


Vo*  a  2 

-j— (Nj 2  +  N2  2 )  +  -9— (N3  2  +  N,2)  +  asan(N1N3  +  N2N4). 


Since  the  quantity  of  interest  is  the  power 
relative  to  the  noise  power  an2 ,  the  equation  is  divided  through 
by  an2.  Letting 


the  signal -to-noise  ratio,  the  relative  power  is 


00  =  %p(N12  +  N22)  +  %(N32  +  N42)  +  vp(N1N3  +  N2N4). 

at 

Since  the  signal  may  be  smeared  over  several 
frequency  bins,  it  is  necessary  to  weight  the  signal -to-noise 
ratio  by  the  fraction  of  the  integration  time  that  the 
target  spends  in  bin  i.  The  expression  for  00  becomes 

00i  =  ky ±p (Nj 2  +  n22)  +  %(N32  +  N42)  /7Ip( n2n3  +  n2n4). 
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It  is  seen  that  the  average  power  relative  to 
noise  power  out  of  the  narrowband  filter  of  an  omnidirectional 
hydrophone  can  be  simulated  by  knowing  the  signal -to -noise 
ratio  and  generating  four  independent  samples  from  a  standard 
normal  distribution.  The  bin  with  the  largest  value  for  00  is 
chosen  by  the  MAX-OR  processor  and  the  center  frequency  of  that 
bin  is  used  as  the  frequency  data  point. 

3 . 3 . 2 . 1  Generation  of  Bearing  Data.  To  arrive  at 

an  estimate  of  bearing,  the  voltage  from  the  omniphone  is 
multiplied  and  time  averaged  with  the  voltage  from  the  X  phone 
and  the  voltage  from  the  Y  phone.  The  bearing  estimate  &  is 
then 


8  -  tan  1  (§£) 

where  OY  and  OX  are  the  averages  mentioned  above. 

The  equations  used  in  the  simulation  for  OX 
shall  be  derived.  The  derivation  of  OY  is  similar.  The  input 
to  the  omniphone  is 

W(t)  =  S(t)  +  n (t) 

as  before.  The  input  to  the  X  phone  is 
V(t)  =  S(t)cos(0)  +  nx(t) 
where  9  is  the  true  target  bearing  and 
nx(t)  -  N(0,  — ) . 


64 


v 


Tracor  Applied  Sciences 


It  can  be  shown  that  n(t)  and  nx(t)  are 
independent  processes.  In  terms  of  orthogonal  components, 

W(t)  =  Sjcosco0t  +  S2sino)0t  +  r^cosojgt  +  n2sinco0t  ; 

V(t)  =  Sj  coso)otcos0  +  S2  sinwotcos0  +  riXicosa>0t  +  nX2sina)0t 

where 

S>  ~  N(0,as2) 

S2  ~  N(0,as2) 
n j  ~  N(0,an2) 
n2  ~  N(0,an2) 
nXi  -  N(0,^~) 
nx2  ~  N(0,^-) . 

W(t)  and  V(t)  are  multipled  and  time  averaged.  It  is  assumed 
that  Si,  S2  ,  m  ,  n2,  nXl  >  and  r'x2  are  slowly  varying.  The 
terms  with  cosoj0tsincu0t  factors  integrate  to  zero.  The 
expression  is  then 

1 

W(t)  V(t)  =  5^[cos2cu0t (S:  2cos6  +  n1nXl  +  +  HjSjCOse) 

+  sin2(u0t (S2 2 cosQ  +  n2nX2  +  S2nX2  +  n2S2cos6)]  dt 
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Integration  yields 


W(t)  V(t)  =  %(S12cos6  +  +  n1Slcos0) 

+  %(S22cos0  +  n2nXz  +  S2nX2  +  n2S2cos9). 

Using  the  transformations  and  the  following 
transformations , 


nxi  = 


_  an 


n 


N. 


nx  =  —  N, 
2  /7  6 


where  Ns ,  Ng  -  N(0,1)  yields 


W(t)  V(t)  =  ^f-d^2  +  N22)cos0  4-  -^-(N3N5  +  N„N6) 

2  2  /7 


+  +  N2N  )  +  +  N2Nlt)cos0 

2/7  2 


Again,  the  quantity  of  interest  is  the  average 
power  relative  to  the  noise  power  a^2  so  the  equation  is  divided 
through  by  a  2 .  Letting 


P  = 


_  aS 


2  ' 
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OX  =  w^t-)-  yCU.  =  £-(N.  2  +  N  2)cos0  +  — — — (H  N  +  N  N.) 

2  2  2  2  n  3  5  4  6 


+  r^p-(N1Ns  +  N2N6)  +  4(NiN3  +  NA>COS0  . 


Once  again,  p  must  be  weighted  by  Y^,  the  fraction 
of  time  spent  in  frequency  bin  i.  The  expression  for  OX  is  now 

OXi  =  +  N2 2 ) cos6  +  — - — (N  N  +  N  N  ) 

2  2  2  /2  35 

Y  _•  P  /  2  /yTp 

+  — (NjN;  +  N2Ns)  +  _i-(N1N3  +  N2Nlt)cose  . 


Y,*  P 

— (N  2  +  N  2)  sin0  +  — - — (N  N  +  NUN, ) 

2  1  2  2  ft  3 

^  Y- o / 2  /yTp 

+  -  ^  i(N1N7  +  N2N8)  +  -~-(N1N3  +  N2N4)sin0  , 

where  N?  ,  N8  -  N(0,1). 

3.3.3  Summary 

The  MAX-OR  processor  chooses  the  frequency  bin 
on  the  basis  of  max{00j}  where  j  ranges  over  the  number  of  bins. 
The  bearing  estimate  is  computed  using  the  values  of  OX  and  OY 
for  that  bin.  The  output  of  the  processor  depicted  in  Figure 
3. 3. 2.1  may  be  modeled  with  the  proper  statistics  by  generating 


Similarly, 


OY.  = 
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eight  samples  of  a  standard  normal  distribution  and  using  the 
above  equations.  The  algorithm  model  of  this  processor  is 
depicted  in  Figure  3.3.3. 

3.3.4  Signal-to-Noise  Ratios  and  Variances 

The  data  used  in  the  simulation  was  generated 
by  an  algorithm  designed  to  generate  data  with  statistical 
properties  closely  resembling  those  of  the  output  of  a  general 
processor  for  DIFAR  buoys.  In  order  to  generate  data  by  this 
scheme,  knowledge  of  the  signal -to -noise  ratio  at  the  input  to 
the  comb  filter  bank  is  required.  The  signal -to-noise  ratio 
was  computed  at  each  time  point  from 

[ S/N] ±  =  S  -  N  -  20  log10  Ri 


where 


[ S/N] ^  is  the  signal -to-noise  ratio  in  dB 
at  time  ^ 

S  is  the  source  level  in  dB 
N  is  the  noise  level  in  dB 
R^_  is  the  range  in  yards  at  t^. 

The  signal  co-noise  ratio  is  input  into  the 
data  generation  algorithm.  This  ratio  is  used  in  the  scheme 
to  generate  the  average  power  in  each  filter  bin  over  the 
integration  time.  At  the  end  of  the  data  generation  algorithm, 
when  a  frequency  bin  has  been  chosen  to  estimate  the  signal 
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FIGURE  3.3.3 


FLOW  CHART  OF  PROCESSOR 
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FIGURE  3.3.3  (Continued) 
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frequency,  an  estimate  of  the  signal -to-noise  ratio  is  made 
using  only  the  average  power  in  the  selected  bin  of  the  comb 
filter  and  the  knowledge  of  the  average  noise  power. 


[ S/N] i  =  10  log10 


Letting  represent  S+N  in  the  selected  bin  and  remembering 

Ill  3.X  ^ 

that  all  power  values  are  relative  to  the  noise  power,  the 
expression  becomes 


S/N  i  =  10  log10 


where  00max  is  measured  in  power  units  (not  dB) .  This  method 
of  estimating  the  signal-to-noise  ratio  was  chosen  to  make  the 
simulation  as  realistic  as  possible. 

The  MLP ,  hybrid,  and  iterated  sequential 
algorithms  require  that  the  covariance  matrix  of  the  data  be 
supplied.  Since  the  data  points  are  independent,  the  covariance 
matrix  is  diagonal  and  all  that  needs  to  be  found  is  an 
estimation  of  the  variance  of  the  population  from  which  each 
data  point  is  obtained.  It  is  clear  that  the  population  and 
the  associated  variances  are  functions  of  the  signal-to-noise 
ratio.  Estimates  of  the  variances  as  a  function  of  signal-to- 
noise  ratio  were  found  by  generating  1,000  samples  for  Af  =  .1  Hz 
at  constant  [S/N]  for  several  different  [S/N]  ratios.  Since  the 
value  of  the  frequency  is  only  resolvable  to  one  bin  width,  the 
frequency  variance  was  given  a  lower  bound  of  variance 

of  a  variate  with  density  function  constant  over  the  interval 
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Af  and  zero  everywhere  else.  These  numbers  are  correct  when 
Af  =  .1  Hz.  For  other  values  of  Af  the  transformation 

aAf  -  (10)  (Af)  (a10) 

is  used  where  Cj 0  is  the  standard  deviation  of  the  frequency 
data  for  a  particular  signal-to-noise  ratio  when  Af  =  .1  Hz. 
The  bearing  variances  do  not  depend  on  Af. 

3 .4  Test  Results  and  Analysis 

3.4.1  Introduction 


This  section  presents  the  results  of  applying 
the  hybrid,  ML.P ,  and  sequential  tracking  algorithms  to  the 
four  scenarios  described  in  Section  3-2.  To  help  understand 
the  effectiveness  of  the  algorithms,  several  plots  and  charts 
have  been  prepared  for  each  scenario.  Figures  3.4.1  through 
3.4.4  contain  the  average  estimated  trajectory  for  each  scenario, 
which  are  found  by  averaging  over  all  Monte  Carlo  runs  the 
position  vectors  for  each  time  point  that  a  state  estimate  is 
generated. 


Figures  3.4.5  through  3.4.8  present  graphs  of 
the  average  distance  error,  over  all  Monte  Carlo  runs,  for  a 
given  scenario.  It  is  computed  by  finding  the  difference 
between  the  estimated  and  true  trajectories  for  each  Monte 
Carlo  run  and  then  averaging  them.  The  dotted  lines  represent 
the  one-sigma  deviation  limits  about  each  average.  Note,  the 
tighter  these  bands  are  around  the  average,  the  smaller  the 
trajectory  error  variation  about  the  average  trajectory  error, 
indicating  that  the  algorithm  has  converged  to  the  true 
trajectory. 
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Figure  3.4.9  gives  the  average  distance  error 
over  the  entire  scenario  (5)  for  all  three  algorithms.  The 
number  to  the  left  of  the  slash  gives  the  actual  average 
distance  error;  the  number  to  the  right  of  the  slash  gives  the 
average  distance  error  as  a  fraction  of  the  hybrid  algorithm's 
distance  error.  Lastly,  Figure  3.4.10  gives  the  average 
execution  time  for  one  Monte  Carlo  run,  computed  by  the 
formula: 


(Total  CPU  Time) / 25 . 

Again,  the  number  to  the  left  of  the  slash  is  the  average 
execution  for  that  algorithm  and  the  number  to  the  right  of 
the  slash  is  that  time  expressed  as  a  fraction  of  the  average 
execution  time  for  the  hybrid. 

3.4.2  Scenario  One 

Figures  3.4.1a-c  and  3.4.5a-c  contain  the 
average  estimated  track  and  average  error  plots  for  scenario 
one.  This  scenario  was  a  simple  straight  line  trajectory 
through  the  center  of  the  buoy  field  and  all  algorithms  followed 
it  well  with  the  following  exceptions  : 

(1)  The  MLP  displays  a  slight  offset  bias. 

This  will  be  evident  in  all  attempts 

of  the  MLP  to  track  straight  line  data. 

(2)  The  iterated  sequential  algorithm  appears 
to  begin  tracking  before  the  start  of  the 
true  track.  This  is  caused  by  the 
inability  of  the  initializer  to  provide 
accurate  state  vector  estimates  with  a 
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small  number  of  data  points,  resulting 
in  severe  underestimation  of  position 
coordinates  early  in  the  scenario. 

This  problem  will  be  evident  in  all 
subsequent  scenarios. 

From  the  error  plots  it  can  be  seen  that  the 
hybrid  has  a  smaller  average  distance  error  than  either  the 
MLP  or  the  sequential,  and  it  also  has  somewhat  tighter  sigma 
lines.  The  error  plots  for  the  sequential  indicate  large 
average  distance  errors  until  about  200  seconds,  reflecting 
the  algorithm's  inability  to  provide  good  position  estimates 
early  in  the  scenario.  Note  that  the  sequential  and  hybrid 
algorithms  begin  generating  estimates  approximately  50  seconds 
into  the  scenario;  the  MLP  cannot  begin  generating  estimates 
until  almost  100  seconds  of  the  scenario  have  passed. 

From  Figure  3.4.9,  the  average  distance  error 
for  the  iterated  sequential  is  about  5.5  times  that  of  the 
hybrid  and  the  average  error  for  the  MLP  is  about  1.8  times 
that  of  the  hybrid.  Lastly,  from  Figure  3.4.10,  a  sequential 
run  took  about  2.8  times  as  long  as  a  hybrid  run  and  an  MLP 
run  took  almost  four  times  as  long.  Thus,  for  Scenario  1,  the 
hybrid  had  smaller  average  tracking  error,  small  tracking  error 
variance,  and  substantially  faster  execution  time.  All  three 
methods  were  within  500  meters  of  the  target  throughout  the 
entire  scenario. 

3.4.3  Scenario  Two 


Figures  3.4.2a-c  contain  the  average  estimated 
track  plots  and  Figures  3.4.6a-c  contain  the  average  distance 
error  plots  for  Scenario  2.  The  hybrid  approximated  the  true 
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trajectory  very  closely,  while  the  MLP  again  exhibited  some  bias 
and  also  had  some  difficulty  with  the  beginning  of  the  turn. 

The  sequential  again  had  estimation  problems  in  the  early  going 
and  then  underestimated  the  turn  somewhat.  The  turn  under¬ 
estimation  occurred  primarily  because  initialization  took  place 
shortly  before  entering  the  turn  (on  the  average)  and  usually 
introduced  some  lag  in  the  state  vector.  This  then  caused  the 
sequential  to  be  somewhat  behind  throughout  the  turn. 

Note  that  the  sigma  limits  for  the  hybrid 
algorithm  are  much  tighter  than  for  the  other  two  algorithms. 

The  troughs  in  the  plot  indicate,  roughly,  those  places  where 
the  batch  portion  of  the  algorithm  had  completed  reinitializa¬ 
tion  of  the  trajectory.  For  the  MLP,  the  large  spike  occurring 
at  approximately  300  seconds  indicates  that  the  turn  has  begun 
and  that  the  MLP  has  lost  the  trajectory  somewhat.  For  the 
sequential  algorithm,  there  are  fairly  substantial  errors 
occurring  at  the  start  and  near  the  end  of  the  scenario.  Both 
sets  of  these  errors  are  due  to  the  inability  of  the  sequential 
starter  to  consistently  provide  accurate  estimates  of  the  state 
vector  using  a  small  number  of  data  points. 

3.4.4  Scenario  Three 


Figures  3.4.3a-c  and  3.4.7a-c  present  the 
average  estimated  track  and  average  distance  error  plots  for 
Scenario  3.  Again,  the  hybrid  estimates  the  track  quite  well 
with  a  stable  error  plot  and  fairly  tight  sigma  limits.  The 
MLP  exhibits  problems  in  two  areas,  recognizing  that  the  turn 
has  begun  and  switching  to  the  appropriate  model,  and,  about 
halfway  through  the  turn,  realizing  that  some  model  adjustment 
is  needed.  Figure  3.4.9  shows  that  the  average  error  for  the 
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MLP  is  about  2.3  times  that  of  the  hybrid.  The  sequential 
again  exhibits  difficulties  for  about  the  first  200  seconds 
of  the  scenario,  causing  it  to  lag  somewhat  throughout  the 
maneuver.  Figure  3.4.10  shows  substantially  faster  execution 
time  and,  again,  all  algorithms  estimated  positions  within 
500  meters  of  the  actual  ones  throughout  the  entire  scenario. 

3.4.5  Scenario  Four 


The  average  estimated  track  and  average  distance 
error  plots  for  Scenario  4  are  contained  in  Figures  3.4.4a-c 
and  3.4.8a-c.  All  three  algorithms  estimate  this  trajectory 
fairly  well,  with  the  hybrid  giving  just  a  bit  lower  average 
error  than  the  MLP.  It  is  interesting  to  note  that  during  the 
maneuver  the  sigma  limits  for  the  hybrid  were  somewhat  larger 
than  for  the  MLP,  but  after  the  maneuver  was  completed  they 
were  substantially  smaller.  All  three  algorithms  required 
more  execution  time  for  this  scenario  than  for  any  other, 
however,  the  sequential  still  took  about  1.8  times  and  the 
MLP  about  3.4  times  as  long  to  run  as  the  hybrid. 
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4.0  CONCLUSIONS  AND  RECOMMENDATIONS 

4.1  Conclusions 


The  main  purpose  of  this  study  was  to  refine 
the  hybrid  algorithm  by  improving  its  numerical  procedures, 
program  organization,  and  switching  rules,  and  then  compare 
its  performance  to  that  of  Tracor's  MLP  algorithm  and  to  that 
of  an  algorithm  which  uses  an  iterated  sequential  filter  to 
initialize.  The  ideas  underlying  the  development  of  all  three 
algorithms  were  described  in  Sections  1.0  and  2.0;  Section  3.0 
presented  the  results  of  a  series  of  tests  designed  to  test 
the  capabilities  of  each  algorithm. 

Based  on  the  results  contained  in  Section  3.4, 
it  is  clear  that  the  hybrid  is  superior  to  both  the  MLP  and  the 
iterated  sequential  in  terms  of  tracking  accuracy  and  execution 
time.  In  every  scenario  examined  the  hybrid  had  not  only  the 
lowest  average  execution  time  (often  by  a  factor  of  three  or 
four),  but  also  the  lowest  average  distance  error.  In  addition, 
the  hybrid  usually  had  tighter  sigma-bounds  about  the  distance 
error  than  either  of  the  other  two. 

Comparisons  between  the  iterated  sequential  and 
the  MLP  are  somewhat  more  difficult  to  make.  For  all  scenarios 
examined,  the  iterated  sequential  executed  considerably  faster 
than  the  MLP.  Also,  once  the  sequential  had  enough  points  to 
provide  correct  initialization,  examination  of  the  distance 
error  plots  indicates  no  significant  differences  between  the 
two  algorithms  in  their  distance  errors  and  distance  error 
variances.  Thus,  it  is  only  early  in  a  given  initialization  or 
reinitialization  phase  that  the  iterated  sequential  does  not 
perform  as  well  as  the  MLP. 
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Part  of  the  reason  that  the  sequential  algorithm 
has  difficulty  initializing  and  reinitializing  is  that  when  the 
extended  Kalman  Filter  measurement  equations  are  solved,  there 
are  no  provisions  for  optimizing  the  resultant  state  vector. 

All  optimization  comes  strictly  as  a  result  of  iteratively 
obtaining  an  initial  guess,  filtering  to  obtain  a  state  estimate, 
and  then  mapping  the  state  vector  back  in  time.  Work  has  begun 
on  the  software  required  to  optimize  the  state  vector  in  terms 
of  the  actual  measurement  equation,  not  the  linearized  approxi¬ 
mation  obtained  when  applying  the  Kalman  gain  equation.  If 
this  procedure  works  as  expected,  both  the  execution  time  and 
the  average  distance  errors  for  the  sequential  will  be  reduced. 

At  present  the  sequential  and  MLP  algorithms  may  be  considered 
roughly  equal  in  performance.  With  the  optimized  initializer, 
the  sequential  will  be  clearly  superior  to  the  MLP. 

4.2  Recommendations 


The  recommendations  arising  from  this  study  fall 
into  two  natural  groups ,  those  dealing  with  the  improvement  or 
modification  of  one  of  the  algorithms  and  those  dealing  with 
analysis  of  the  capabilities  and  robustness  of  each  algorithm. 
Recommendations  for  algorithm  improvement  or  modification 
include : 

(1)  Dealing  with  Outliers.  At  present, 

neither  the  hybrid  nor  the  sequential 
has  an  effective  method  for  dealing 
with  bad  data  points  or  a  bad  data 
stream  from  one  particular  sensor. 

There  are  several  schemes  for  dealing 
with  outliers  (such  as  the  procedure 
used  in  the  MLP)  and  they  should  be 
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investigated  for  use  in  the  hybrid  and 
sequential  algorithms. 

(2)  Optimize  Sequential  Initializer.  As 
noted  in  Section  4.1,  software  is  being 
developed  which  will  allow  the  sequential 
starter  to  optimize  a  given  state  vector 
estimate  with  respect  to  the  true  measure 
ment  models  and  not  some  linear  approxi¬ 
mation.  This  should  be  completed  and  the 
new  optimal  initializer  installed. 

(3)  Add  Higher  Order  Terms .  For  the  hybrid 
algorithm,  the  optimization  procedure 
used  in  the  batch  initializer  was 
equivalent,  in  result,  to  using  higher 
order  terms  in  the  measurement  model 
approximations.  However,  in  the 
sequential  portion  of  the  algorithm 
these  optimization  procedures  were  not 
implemented  in  order  to  minimize 
execution  time.  By  using  higher  order 
terms  in  the  measurement  models  employed 
by  the  sequential  filter,  it  may  be 
possible  to  increase  tracking  accuracy 
without  significantly  raising  scenario 
execution  time. 

Recommendations  for  determining  and  comparing 
the  robustness  of  each  algorithm  include: 
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(1)  Tests  on  Real  Data.  Using  the  non- 
gaussian,  simulated  data  of  this  study, 
all  three  tracking  algorithms  did  well. 
However,  the  true  test  of  any  algorithm's 
capabilities  is  its  performance  on  actual 
sea  data.  Once  the  modifications  proposed 
above  have  been  completed,  all  three 
algorithms  should  be  tested  on  real  data 
taken  from  several  scenarios. 

(2)  Effects  of  Data  Quality  and  Data  Rate. 
During  testing  of  the  hybrid,  scenarios 
with  different  data  rates  and  qualities 
were  produced  and  tested.  Indications 
were  that  data  rate  was  more  a  factor 
in  determining  good  tracking  accuracy 
than  was  data  quality.  Using  analysis 
of  variance/response  surface  techniques, 
it  may  be  possible  to  determine  the 
relative  importance  of  data  rate,  data 
quality,  and  buoy  number  and  to  identify 
certain  optimal  conditions. 

(3)  Effect  of  Buoy  Drift.  All  scenarios 
analyzed  in  this  study  assumed  constant 
buoy  positions  throughout  the  scenario. 

Of  course,  in  actual  practice  this  is 
not  the  case  and  the  best  algorithm  is 
that  one  which  is  most  efficient  in  the 
face  of  buoy  position  uncertainty.  The 
necessary  software  is  in  place  in  the 
hybrid  and  sequential  to  generate  buoy 
position  estimates,  however,  there 
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would  be  a  certain  amount  of  programming 
involved  in  generating  the  simulated 
measurements . 

(4)  Examine  Other  Data  Types.  For  this  study 
the  only  data  types  used  for  each  scenario 
were  frequency  and  bearing.  However, 
there  are  other  data  measurements  which 
can  be  used  for  tracking,  such  as  range, 
time  difference  of  arrival,  Doppler  ratio, 
and  Doppler  difference. 
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5.0  MULTIPLE  TARGET  PROBLEM 


The  second  major  purpose  of  this  study  was  to 
outline  a  procedure  for  attacking  the  multiple  target  estimation 
problem.  There  are  two  facets  to  the  problem: 

(1)  Having  the  ability  to  track  one  particular 
target  in  the  midst  of  others. 

(2)  Having  the  ability  to  identify  and  track 
several  different  targets  concurrently. 

These  may  be  considered  to  be  the  surveillance  and  tracking 
aspects  of  the  multiple  target  problem.  It  is  hoped  that  a 
single  algorithm  can  be  developed  to  handle  both  tracking  and 
surveillance . 


At  present  there  are  several  theoretical  papers 
(References  8-15)  dealing  with  various  aspects  of  the  problem. 
The  majority  of  them,  however,  do  not  consider  in  depth  the 
problem  of  initial  sorting  and  classification,  concentrating 
instead  on  schemes  for  placing  an  observation  in  the  correct 
track,  given  that  a  certain  number  of  tracks  exist.  There  are 
several  related  approaches  which  basically  develop  the  theory 
needed  to  set  up  probability  "gates"  around  the  predicted 
measurement  for  a  given  track.  These  gates  are  equiprobability 
contours  chosen  such  that  the  likelihood  of  any  estimate  within 
the  gate  being  correct  is  above  a  certain  threshold.  Other 
methods  rely  on  a  posteriori  probability  analyses  of  the  likeli¬ 
hood  of  a  given  measurement  belonging  to  a  certain  track,  with 
appropriate  decision  rules  for  placing  the  observation  with  a 
particular  track.  Many  papers  contain  algorithms  for  multiple 
target  tracking  which  require  either  prohibitive  amounts  of 
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computer  core  or  prohibitive  amounts  of  computer  time, 
necessitating  the  use  of  some  approximating,  suboptimal 
procedure.  Currently,  no  paper  has  yet  been  found  which 
applies  a  given  method  to  a  real-world  problem.  All  examples 
given  are  fairly  simple  simulations. 

It  is  proposed  to  approach  this  problem  by 
using  a  combination  of  the  above  methodologies  and  a  statistical 
technique  called  cluster  analysis.  Cluster  analysis  is 
essentially  the  branch  of  statistics  which  deals  with  methods 
for  grouping  large  numbers  of  objects  into  smaller,  mutually 
exclusive  subgroups  containing  members  as  much  alike  as  possible. 
Many  clustering  methods  are  extremely  effective  when  applied  to 
the  appropriate  data  set.  Clustering  procedures  may  all  be 
thought  of  as  containing  three  basic  elements: 

(1)  A  measure  of  similarity  (or  dissimilarity) 
to  apply  to  members  of  the  population. 
Examples  include  Euclidean  distance, 
weighted  Euclidean  distance,  and  the 
value  of  some  scoring  function. 

(2)  Some  optimizing  or  cluster  defining 
criterion.  Examples  include: 

(a)  Minimizing  the  variance  of  each 
cluster  about  its  centroid. 

(b)  Minimizing  the  maximum  distance 
between  any  two  cluster  points  for 
all  possible  clusters. 
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(c)  Minimizing  the  average  distance 
between  all  cluster  points  over 
all  possible  clusters. 

(3)  An  algorithm  or  method  for  finding  the 
optimum  based  on  (1)  and  (2)  above.  In 
general  these  algorithms  fall  into  two 
classes : 

(a)  Divisive  or  agglomerative  partitioning 
procedures,  such  as  sorting,  switching, 
joining,  splitting,  or  adding. 

(b)  Application  of  some  standard 
mathematical  formulation,  such  as 
integer  programming ,  dynamic 
programming,  or  graph  coloring. 

Figure  5.1  presents  a  logical  flowchart  of  the 
proposed  multiple  target  algorithm.  The  algorithm  would  contain 
a  clustering  procedure  which  would  perform  two  functions: 

(1)  Initially  break  a  set  of  data  points 
into  clusters  corresponding  to  a  set 
of  trajectories. 

(2)  Process  each  data  point  after  initializa¬ 
tion  to  determine  which  cluster  (trajectory) 
it  should  be  associated  with. 

Once  a  point  had  been  associated  with  a  given  cluster,  a  new 
state  vector  estimate  would  be  generated  using  either  the  hybrid 
or  iterated  sequential  algorithm.  Note  that  the  speed  of  the 
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particular  tracking  algorithm  becomes  crucial  at  this  time, 
particularly  as  the  number  of  targets  increases.  The  resulting 
state  vector  would  then  be  analyzed  for  appropriateness  using 
another  clustering  algorithm,  some  type  of  a  posteriori 
probability  analysis,  or  one  of  the  "gating"  procedures 
discussed  in  the  literature.  If  an  observation  passes  this 
test  it  becomes  permanently  associated  with  the  current  cluster 
and  a  new  observation  is  obtained.  If  an  observation  fails  the 
appropriateness  test,  it  is  placed  in  the  next  most  likely 
cluster  (that  cluster  having  the  next  highest  value  of  the 
optimizing  criteria)  and  a  new  state  vector  estimate  is  derived 
based  on  the  trajectory  associated  with  members  of  this  cluster. 
If  post -estimation  tests  rule  an  observation  out  of  all  its 
most  likely  clusters,  a  new  cluster  is  formed  with  the  given 
observation  as  its  sole  member.  This  cluster  will  be  considered 
the  beginning  of  a  new  trajectory.  Note  that  the  state  vector 
estimation  algorithm  is  being  used  in  an  interactive  manner. 

It  is  hoped  that  by  having  both  pre -estimation  and  post¬ 
estimation  evaluations  of  each  observation  the  number  of 
misclassif ications  can  be  kept  to  a  minimum. 

Because  of  the  difficulty  of  the  problem  it  is 
felt  that  all  data  processing  procedures  should  be  made  general 
enough  to  handle  N  targets,  but  that,  initially,  all  efforts 
should  concentrate  on  being  able  to  handle  one  to  three  targets 
in  terms  of  tracking  or  surveillance.  When  this  stage  is 
reached,  an  analysis  of  algorithm  performance  and  a  study  of 
the  additional  effort  needed  to  handle  N  targets  should  be  made. 

At  this  stage,  then,  identifiable  tasks  in  a 
multi-target  project  include: 
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(1)  Generation  of  simulated  multi-target 
data. 

(2)  Selection  and  testing  of  an  appropriate 
clustering  algorithm. 

(3)  Selection  and  testing  of  an  appropriate 
post -estimation  evaluation  procedure. 

(4)  Development  of  software  to  support 
and  implement  the  multiple  target 
algorithm. 

(5)  Modification  of  hybrid  and/or  iterated 
sequential  tracking  algorithm  to 
interface  with  multi-target  program. 

This  would  involve  some  storage 
alterations  and  possible  some  new 
subroutines  to  handle  efficiently  the 
data  structures  that  the  multi-target 
program  will  require. 

Certain  parts  of  the  above  schedule  have  already 
been  addressed,  either  in  this  project  or  other  Tracor  projects 

(1)  For  this  project  the  capability  to 

generate  simulated  frequency  and  bearing 
data  for  one  target  was  developed.  It 
is  felt  that  two  targets  can  be  simulated 
by  generating  two  sets  of  data  and  merging 
the  results.  A  small  amount  of  programming 
will  be  required  to  develop  the  software 
required  for  the  merge . 


Tracor  Applied  Sciences 


(2)  Tracor  is  now  in  the  process  of 
acquiring  two  general  purpose  clustering 
routines,  CLUSTAR  and  CLUSTID  (Reference  15), 
which,  according  to  its  authors,  will  allow 
the  examination  of  about  ”75  percent  of  the 
published  uses  of  cluster  analysis". 
Additionally,  Tracor  is  in  the  process  of 
obtaining  an  algorithm  developed  by 

R.  F.  Ling  (References  16,  17)  which  is 
specifically  designed  to  produce  long, 
string-like  clusters  very  similar  to  the 
graphs  of  time,  frequency,  and  bearing 
coordinates  of  a  moving  target.  Ling 
used  it  very  successfully  to  cluster 
the  60  brightest  stars  in  the  sky  into 
their  respective  constellations  using 
celestial  coordinates. 

(3)  A  literature  search  is  continuing  for 
additional  papers  specifically  addressing 
the  multi-target  problem. 
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